#ifndef _FEMSPACE_UNSTRUCT_H_
#define _FEMSPACE_UNSTRUCT_H_
#include "FEMSpace.h"
#include <cassert>

/**
 * @brief 返回map，存放边界上自由度的全局标号，和边界函数(仅Dirichlet)在对应自由度上的函数值。
 *
 * @param _boundary_function 边界函数
 * @param _boundary_dofs 边界上的自由度
 */
template <unsigned int DIM>
std::map<int, double> interpolateBoundaryValues_Unstruct(Mesh<DIM> *_mesh, vector<function<double(double *)>> &_boundary_function, std::set<Dofs<DIM>> &_boundary_dofs)
{
    EasyMesh *easy_mesh = dynamic_cast<EasyMesh *>(_mesh);
    if (easy_mesh == nullptr)
    {
        cout << "casting failed" << endl;
        exit(EXIT_FAILURE);
    }

    std::map<int, double> boundary_values;
    for (auto &dofs : _boundary_dofs)
    {
        int index = dofs.getGlobalIndex();
        // EasyMesh中根据自由度index取出boundarymark
        int boundarymark = easy_mesh->getBoundaryMark(index);
        double p[2];
        for (int i = 0; i < 2; ++i)
        {
            p[i] = dofs[i];
        }
        double value = _boundary_function[boundarymark - 1](p);
        boundary_values.insert({index, value});
    }
    return boundary_values;
};

/**
 * @brief 对刚度矩阵和右端项应用Dirichlet条件
 *
 * @param system_matrix 刚度矩阵
 * @param system_rhs 右端项
 */

class Possion_2D_Unstruct : public FEMSpace<2>
{
public:
    Possion_2D_Unstruct();
    Possion_2D_Unstruct(Mesh<2> *, TemplateElement<2> *, Equation<2> *, Domain<2> *);
    Possion_2D_Unstruct(const Possion_2D_Unstruct &);
    ~Possion_2D_Unstruct(){/*delete _mesh;delete _element;delete _equation;*/};
    // virtual void make_grid();
    virtual void setup_system();
    virtual void assemble_system();
    virtual void dealwithBoundaryCondition();
    virtual void solve();

    double LocalStiffMatrix(int i, int j);
    Solution getSolution() const;
    SpMat A() const;
    RHS Rhs() const;

protected:
    SpMat system_matrix;
    RHS system_rhs;
    Solution solution;
    Domain<2> *domain;
    Mesh<2> *mesh;
    TemplateElement<2> *element;
    Equation<2> *equation;
};

double Possion_2D_Unstruct::LocalStiffMatrix(int i, int j)
{
    double a = 0;
    for (int k = 0; k < element->n_GaussPnt(); k++)
    {
        double xi = element->GaussionPoint(k)[0];
        double eta = element->GaussionPoint(k)[1];
        double b = 0;
        for (int l = 0; l < element->gradient(xi, eta, i).size(); l++)
            b += element->gradient(xi, eta, i)[l] * element->gradient(xi, eta, j)[l];
        a = a + element->det_Jacobi(xi, eta) * element->GaussionWeight(k) * b;
    }
    return a;
}

Solution Possion_2D_Unstruct::getSolution() const
{
    return solution;
}

Possion_2D_Unstruct::Possion_2D_Unstruct()
{
    // TODO:如果不是P1单元，这样写目前是不对的
    int num = mesh->getTotalNumDofs();
    system_matrix = SpMat(num, num);
    system_rhs = RHS(num);
}
Possion_2D_Unstruct::Possion_2D_Unstruct(Mesh<2> *_mesh, TemplateElement<2> *_element, Equation<2> *_equation, Domain<2> *_domain)
{
    domain = _domain;
    mesh = _mesh;
    element = _element;
    equation = _equation;
    // TODO:如果不是P1单元，这样写目前是不对的
    int num = mesh->getTotalNumDofs();
    system_matrix = SpMat(num, num);
    system_rhs = RHS(num);
}

Possion_2D_Unstruct::Possion_2D_Unstruct(const Possion_2D_Unstruct &p)
{
    domain = p.domain;
    mesh = p.mesh;
    element = p.element;
    equation = p.equation;
    system_matrix = p.system_matrix;
    system_rhs = p.system_rhs;
}

void Possion_2D_Unstruct::setup_system()
{
    /**
     * set the sizes of system_matrix, right hand side vector ans the solution vector;
     */
    int num = mesh->getTotalNumDofs();
    std::cout << "Number of degrees of freedom: " << num << std::endl;
    // TODO:如果不是P1单元，这样写目前是不对的
    system_matrix = SpMat(num, num);
    system_matrix.setZero();
    system_rhs = RHS::Zero(num);
    solution = Solution::Zero(num);
}

void Possion_2D_Unstruct::assemble_system()
{
    /**
     * Get the right_hand_side function and boundary function
     */
    auto right_hand_side = equation->RightHandsTermFunction();

    /**
     * for each Grid, we construct a cell_matirx and cell_rhs.
     */
    const unsigned int dofs_per_cell = element->n_Dofs();
    vector<Eigen::Triplet<double>> TriList(mesh->getNumGrid() * dofs_per_cell * dofs_per_cell);
    vector<Eigen::Triplet<double>>::iterator it = TriList.begin();
    Eigen::MatrixXd cell_matrix(dofs_per_cell, dofs_per_cell);
    RHS cell_rhs(dofs_per_cell);

    /*assemble the global matrix and vector cell-by-cell.*/
    for (int index = 0; index < mesh->getNumGrid(); ++index)
    {
        // bind the TemplateElement with current cell.
        Grid<2> cell = mesh->getGrid(index);
        element->reinit(cell, *mesh, index);

        cell_matrix = Eigen::MatrixXd::Zero(dofs_per_cell, dofs_per_cell);
        cell_rhs = RHS::Zero(dofs_per_cell);

        for (int i = 0; i < dofs_per_cell; i++)
        {
            for (int j = 0; j < dofs_per_cell; j++)
            {
                *it = Eigen::Triplet<double>(element->getGlobalIndex(i), element->getGlobalIndex(j), LocalStiffMatrix(i, j));
                it++;
            }
        }
        for (int q_index = 0; q_index < element->n_GaussPnt(); ++q_index)
        {
            Point<2> pnt = element->GaussionPoint(q_index);
            double xi = pnt[0], eta = pnt[1];
            double x = element->Global_x(xi, eta);
            double y = element->Global_y(xi, eta);
            double p[2] = {x, y};
            for (int i = 0; i < dofs_per_cell; ++i)
            {
                cell_rhs(i) += (element->phi(xi, eta, i) *
                                right_hand_side(p) *
                                element->det_Jacobi(xi, eta) *
                                element->GaussionWeight(q_index));
            }
        }
        for (int i = 0; i < dofs_per_cell; ++i)
            system_rhs[element->getGlobalIndex(i)] += cell_rhs(i);
    }
    system_matrix.setFromTriplets(TriList.begin(), TriList.end());
    system_matrix.makeCompressed();
    // cout << system_matrix << endl;
}

void Possion_2D_Unstruct::dealwithBoundaryCondition()
{
    vector<function<double(double *)>> boundary_function = equation->BoundaryConditionFunction();
    cout << "boundary function的size为:" << boundary_function.size() << endl;
    std::set<Dofs<2>> boundary_dofs;
    boundary_dofs = mesh->getBoundaryDofs();
    for (auto dof = boundary_dofs.begin(); dof != boundary_dofs.end(); ++dof)
    {
        cout << *dof << "\n";
    }
    cout << endl;
    // TODO:采用新函数计算边界自由度上的函数值。

    cout << "interpolateBoundaryValue\n"
         << endl;
    // std::map<int, double> boundary_values;
    // boundary_values = interpolateBoundaryValues_Unstruct(mesh, boundary_function, boundary_dofs);
    //  EasyMesh* easy_mesh = dynamic_cast<EasyMesh *>(mesh);
    //  if (easy_mesh == nullptr)
    //  {
    //    cout << "casting failed" << endl;
    //    exit(EXIT_FAILURE);
    //  }
    EasyMesh *easy_mesh = dynamic_cast<EasyMesh *>(mesh);
    // int bndmark = easy_mesh->getBoundaryMark(0);
    std::map<int, double> boundary_values;
    for (auto &dofs : boundary_dofs)
    {
        double p[2];
        for (int i = 0; i < 2; ++i)
        {
            p[i] = dofs[i];
        }
        cout << "自由度信息:" << dofs << endl;
        cout << "getGlobalindex\n";
        int index = dofs.getGlobalIndex();
        // EasyMesh中根据自由度index取出boundarymark
        cout << "getBoundaryMark:\n";
        int boundarymark = easy_mesh->getBoundaryMark(index);
        assert(boundarymark > 0);
        cout << "get value\n";
        // BUG: why???
        cout << "boundary mark is" << boundarymark << endl;
        double value = boundary_function[boundarymark - 1](p);
        cout << "insert into boundary_values\n";
        boundary_values.insert({index, value});
    }

    applyBoundaryValues(boundary_values, system_matrix, system_rhs);
}

void Possion_2D_Unstruct::solve()
{
    /*<  CG solve */
    Eigen::ConjugateGradient<SpMat> Solver_sparse;
    Solver_sparse.setTolerance(1e-40);
    Solver_sparse.compute(system_matrix);
    solution = Solver_sparse.solve(system_rhs);
    /*<  AMG solve */
    // AMG_Multigrid AmgSolve(system_matrix, system_rhs, 1e-12, GS, 2.0 / 3, 0.5);
    // AmgSolve.Solve();
    // solution = AmgSolve.GetResult();
}

SpMat Possion_2D_Unstruct::A() const
{
    return system_matrix;
}
RHS Possion_2D_Unstruct::Rhs() const
{
    return system_rhs;
}

#else
#endif
